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Abstract 

We consider a class of Boussinesq systems of Bona- Smith type in two space dimensions approxi- 
mating surface wave flows modelled by the three-dimensional Euler equations. We show that various 
initial-boundary-value problems for these systems, posed on a bounded plane domain are well posed 
locally in time. In the case of reflective boundary conditions, the systems are discretized by a modi- 
fied Galerkin method which is proved to converge in L 2 at an optimal rate. Numerical experiments 
are presented with the aim of simulating two-dimensional surface waves in complex plane domains 
with a variety of initial and boundary conditions, and comparing numerical solutions of Bona-Smith 
systems with analogous solutions of the BBM-BBM system. 

1 Introduction 

In this paper we will study the Boussinesq system 

m + V • v + V • nv - bAri t = 0, 

(1.1) 

v t + V?7 + ±V|v| 2 + cVAr? - 6Av t = 0, 

where b > and c < are constants. This system is the two-dimensional version of a system in one space 
variable originally derived and analyzed by Bona and Smith, |BS] . It belongs to a family of Boussinesq 
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systems, |BUSIj . |BCLj . that approximate the three-dimensional Euler equations for the irrotational free 
surface flow of an ideal fluid over a horizontal bottom. In (jl.ip the independent variables x = [x, y) 
and t represent spatial position and elapsed time, respectively. The dependent variables r\ = rj(x, t) and 
v = v(x, t) = (u(x, t), u(x, t)), represent quantities proportional, respectively, to the deviation of the free 
surface from its level of rest, and to the horizontal velocity of the fluid particles at some height above the 
bottom. In (jl.ip the variables are nondimensional and unsealed and the horizontal velocity v is evaluated 
at a height z = — 1 + 9(1 + ?7(x, t)), for some 6 6 [-y/2/3, 1]. (In these variables the bottom lies at height 
z = — 1.) In terms of the parameter 9, the constants in (jl.ip are given by the formulas b = (39 2 — l)/6 
and c = (2 — 3# 2 )/3. (The system originally analyzed by Bona and Smith in BSJ corresponds to 6 2 = 1.) 
The value 9 2 = 2/3 (i.e. c = 0) yields the BBM-BBM system, [BC], [DMSIIj . [Chj . 

The system is derived from the Euler equations, BCSIJ, [BCL , under a long wavelength, small 
amplitude assumption. Specifically, one assumes that e := A/ho "C 1, X/ho ^> 1 with the Stokes number 
S := AX 2 /Hq being of 0(1). (Here A is the maximum amplitude of the surface waves measured over 
the level of rest, z = — ho is the (constant) depth of the bottom, and A is a typical wavelength.) If one 
takes 5 = 1, then, in nondimensional, scaled variables, appropriate asymptotic expansions in the Euler 
equations yield equations of the form 

rj t + V • v + e(V • r/v - bArj t ) = 0(e 2 ), ^ 
v f + Vr; + e (| V|v| 2 + cVA?? - 6Av t ) = 0{e 2 ), 

from which (jl.ip follows by unsealing to remove s, and replacing the right-hand side by zero. 

The Cauchy problem for (jl.ip in the case of one spatial variable has been proved to be globally well 
posed for 2/3 < 2 < 1 in appropriate classical and Sobolev space pairs, [BSJ . [BCSI] . The analogous 
problem for 9 2 = 2/3 is locally well posed, jBCSIj . In UDMSIJ we considered a more general class of systems 
in the two-dimensional case and proved that the corresponding Cauchy problem is locally well posed in 
appropriate pairs of Sobolev spaces. Initial-boundary-value problems (ibvp's) for (jl.ip for 2/3 < 9 2 < 1 
on a finite interval in one space variable were analyzed in [ADMIj . and in |BC| for 9 2 = 2/3. It was proved 
in [ADMIj that the ibvp with Dirichlet boundary conditions, wherein 77 and u are given functions of t 
at the endpoints of the interval, is locally well posed. The corresponding ibvp with reflection boundary 
conditions (rj x = u = at both endpoints of the interval) was shown in }ADMI| to be globally well posed; 
so is also the periodic ivp. In DMSII we analyzed three ibvp's for the BBM-BBM system on a smooth 
plane domain O, corresponding to homogeneous Dirichlet boundary conditions for rj and v on <9f2, to 
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homogeneous Neumann boundary conditions for r\ and v on dQ, and to the (normal) reflective boundary 
conditions f? = 0, and v = on dfl, where n is the normal direction to the boundary. We showed that 
these ibvp's are well posed locally in time in the appropriate sense. 

In Section 2 of the paper at hand we consider the Bona-Smith system (jl.lj) and pose it as an ibvp on a 
plane domain Q under a variety of homogeneous boundary conditions on <9f2, including e.g. homogeneous 
Dirichlet b.c.'s for -q and v and reflective b.c.'s. We prove that the corresponding ibvp's are well posed, 
locally in time. 

Turning now to the numerical solution of ibvp's for systems of the type (jl.ip by Galerkin-finitc 
element methods, we note first that it is quite straightforward to construct and analyze such schemes 
for the BBM-BBM system. For example, in [DMSI| we proved optimal-order L 2 -error estimates for the 
standared Galerkin semidiscretization of the BBM-BBM system with homogeneous boundary conditions 
on a smooth domain with a general triangulation. When c > 0, i.e. in the case of the proper Bona-Smith 
systems, the presence of the term VA77 complicates issues. In |DMSIj we analyzed the standard Galerkin 
semidiscretization with bicubic splines for this class of systems posed on rectangles with homogeneous 
Dirichlet boundary conditions, and proved optimal-order i? 2 -error estimates for the approximation of r\. 
(Experimental evidence indicates that the L 2 -errors for the approximation of 77, u and v with this scheme 
have suboptimal - 0(h 3 ) - rate of convergence. In the one-dimensional case one may derive optimal-order 
estimates for the approximations of 77, u, v in W 1,00 x L°° x L°°, cf. [ADMIIJ.) 

In Section 3 of the present paper we consider the Bona-Smith systems with c posed on convex 
smooth planar domains with reflective boundary conditions. The systems are discretized on an arbitrary 
triangulation of the domain using a modified Galrkin method, wherein the Laplacian in the VA77 terms 
in (jl.ip is discretized weakly by an appropriate discrete Laplacian operator. This enables us to prove 
optimal-order L 2 - and H - error estimates on finite element subspaces of continuous, piecewise polynomial 
functions that include the case of piecewise linear functions utilized in most applications. The systems of 
ordinary differential equations representing the semidiscretizations of the Bona-Smith systems are shown 
to be non-stiff. One may thus use any explicit method for their temporal discretization. 

We close the paper by showing the results of a series of numerical experiments of simulations of 
surface waves in complex domains, aimed at comparing the numerical solution of a Bona-Smith system 
with analogous results obtained by solving the BBM-BBM system. 
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2 Well-posedness of ibvp's for the Bona-Smith system 

Let n be a bounded plane domain with smooth boundary (or a convex polygon). We consider the 
Bona-Smith system 

rjt + V • (v + ryv) - bArj t = 0, (2.1a) 

v t + V (r) + i |v| 2 + cAri) - bAv t = 0, (2.1b) 

for (x, t) G f2 x R + , with 6 > 0, c < 0, with initial conditions 

77(x,0) =»7o(x), v(x,0)=v (x) xGfi. (2.2) 

We now describe the class of boundary conditions on dfl under which the problem (2.1)- (|2.2|) will 
be solved. Let H s = H s (il), set, denote the L 2 -based, real- valued, Sobolev classes on fl and Hq the 
subspace of H 1 whose elements have zero trace on dfl. In the sequel, we shall denote by || • || and (■, ■) 
the norm and inner product, respectively, of L 2 = L 2 (fl), by || • || s the norm of H s , and by || • ||oo the 
norm of L°° = L°°(f2). The boundary conditions on r\ will be of the form 

Bi] = 0, xeffl, tel + , (2.3a) 

where the linear operator B and the domain f2 will be assumed to be such that the boundary-value 
problem 

—bAw + w = /, in f2, 

< 

Bw = 0, in dn, 

i. 

has for each / in L 2 a unique solution w E H 2 for which Bw\qq, is well defined. (We will also assume 
that an H 1 solution w of the problem is defined whenever / G H^ 1 ). The boundary conditions on v will 
be of homogeneous Dirichlet type, i.e. 

v = 0, x G dn, te E+. (2.3b) 

We note that examples of suitable boundary conditions of the type (|2.3a|) include, among other, homo- 
geneous Dirichlet (rj = 0), Neumann (^ = 0) or Robin (af^ + f3rj — 0) boundary conditions on the 
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boundary dfl if fl is a bounded plane domain or a convex polygon, boundary conditions of the form 
§^|ri = 0; v\t2 — f° r a multiply connected domain, for example such as the one shown in Figure [2~Tl 
et al. 




Figure 2.1: Plane domain f2 with 



0, ry|r 2 =0. 



In the sequel we let X := 77 2 (f2) n {w : fiw = on dfl} and Hq = (Hq) 2 etc. The main result of this 
section is 

Theorem 2.1 Given r/o G X, Vo G Aqj ^ ere exists T > and a unique solution (rj,v) G C([0, T];X) n 
C([0,T];HJ) o/tfie ifeup ^.Jj, i fEl) . (2.3). Moreover, for each integer (§^, G C([0,T];X)n 
C([0,i];H o ). 



Proof: Write (|2~Taj) and ([2Tb| as 



?7 t + (/-6A)- 1 V-(v + nv) = 0, 
v t + (7 - fcAJ^V U + i|v| 2 + cArA = 0, 



(2.4a) 
(2.4b) 



for x G CI, t > 0. In (|2~4^|) (7 - &A)" 1 denotes the inverse of the operator 7 — bA with domain X , while 
in (|2.4b[) (7 — 6A) _1 represents the inverse of 7 — b A with domain H 2 n Hq. Let F be the vector field on 
X x HJ defined by 



Ffa, v) := ( (7 - 6A)- X V • (v + rpr), (7 - &A) _1 V(»7 + i|v| 2 + cAn) 



F is well defined on X x Hq, since, by the Sobolev imbedding theorem, ?yv G H 1 and |v| 2 G L 2 . Hence, 
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(I - 6A)- X V • (v + jjv) e X and (/ - bA^Vfa + i|v| 2 + cAn) G Hj. Moreover, F is C 1 on X x Hj, 
with derivative -F' (77*, v*) given by 

F'(t]*,v*)(t], v) = ((/ - V • (v + rjv* + 77V), (I - bAy 1 V(r? + v* • v + cAr?)) . 

The continuity of _F' follows from the Sobolev imbedding theorem and the regularity properties of the 
operators (I — 6A) _1 , considered as inverses of (/ — bA) on X for the first component, and of I — bA on 
H 2 fl Hj for the second. 

By the standard theory of ordinary differential equations in Banach spaces, we conclude therefore that 
there exists a unique maximal solution (77, v) G C 1 ([0,T];X) x C^fO, T]; Hj) of (|2.4p for some T > 0, 
with rj\ t= o = 770 and v\t=o = vq. The first conclusion of the theorem follows. The assertion on ^pr 
follows by differentiating (|2.4[) with respect to t k — 1 times. □ 

Remark 2.1 It is not hard to see, using the energy method on the nondimensional but scaled system 
(|1.2p (with right-hand side replaced by zero), that at least in the cases of homogeneous Dirichlet boundary 
conditions for 77 and v and reflective boundary conditions |^ = 0, v = on dfl, the maximum existence 
time T e is independent of e. 

Remark 2.2 If the domain is smooth, one gets smooth solutions from smooth data. Namely, let 
ken, k > 3, and assume that 7/ G X n H k , v £ Hj n H* -1 , then (77, v) G C([0,T];X n H k ) n 
C([0, T]; Hj fl H fc_1 ). This follows directly from the elliptic regularity estimates on / — bA (with the ad 
hoc boundary conditions). One is thus reduced to an ODE in X Pi H k x Hj fl H fc_1 . 

Remark 2.3 Consider the ibvp (|2.ip - (|2.2l) for the Bona-Smith system with (normal) reflective bound- 
ary conditions 

^ = 0, v = 0, for x G dil, t G R+. (2.5) 
on 

It is known that in one dimension, cf. ADMI , the Hamiltonian of the ibvp (|2.ip . (|2.2p . (|2.5p . 

= ff(77, v) = 1 J [ V 2 + (1 + 77)|v| 2 - c|V7?| 2 ] dx, (2.6) 

is conserved. (Indeed this implies global existence-uniqueness of the solution of this ibvp in ID provided 
that 770(2;) + 1 > and H(rjo, uq) is suitably restricted.) In the two-dimensional case H is not conserved 
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in general. To see this, write the system (|2.ip for (x, i)6fix R + as 



■qt + V • P = 0, 



(2.7a) 



vj + VQ + WxWj^O 



(2.7b) 




Jn 

Using now the reflective boundary conditions in (|2.8|) and integrating by parts we see that in the maximal 
temporal interval of existence of a solution of the ibvp (|2.ip . (12. 2p . (12. 5p . 



Hence, a simple sufficient condition for the conservation of the Hamiltonian is irrotationality of the flow. 
In one space dimension, the flow is trivially irrotational. In 2D taking the curl of (|2.1ap we see that 



Now, if the flow is, for example, irrotational at t = 0, the reflective boundary conditions do not allow 
us to conclude from (|2.10p that w = for t > 0. However, in the case of the Cauchy problem or the 
ibvp with periodic boundary conditions on r\ and v on a rectangle (wherein (|2.8p holds as well), lu(0) = 
implies by ([2~TUf that u> = for t > 0. (This was also noticed for the BBM-BBM system in [CI] .) In 
these cases, it follows by (|2.9p that the Hamiltonian is invariant. Note however that conservation of H 
does not imply global well-posedness in 2D. 

3 A modified Galerkin method 

We turn now to the numerical solution of the ibvp (|2.ip . (|2.2p . (|2.5p by Galerkin- finite clement methods. 
The usual Galerkin method for Bona-Smith systems that was analyzed in DMSI requires finite element 
subspaces consisting of C 2 functions. Hence, it is not very useful in practice, as it cannot be applied to 
arbitrary plane domains and triangulations. In this section a modified Galerkin method for the numerical 




(2.9) 



d t (u - bAuj) = 0, <> 0. 



(2.10) 
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solution of the ibvp with (normal) reflection boundary conditions is analyzed. The method is more 
versatile, being applicable to general triangulations and domains (even with piecewise linear continuous 
functions) and is shown to have error estimates with optimal convergence rates in L 2 and H . 

For ease of reference we rewrite here the ibvp that we will approximate. We seek 77 and v, defined for 
(x, t) E il x [0, T] and satisfying for constants b > and c < 

m + V • v + V • 77V - bAm = 0, 

(x,t) 6flx [0,T], 

vt + V(77 + i|v| 2 +cA7 ? )-6Av t =0, 

(11) 

r)(x, 0) = ?7o(x), v(x, 0) = v (x), x e O, 
ga(x,t) = 0, v(x,i) = 0, (x,t) eaflx [0,T]. 

We assume that 17 is a convex domain with smooth enough boundary dfl and that the ibvp (jS]) has a 
unique solution (77, v) which is smooth enough for the purposes of its numerical approximation. In the 
sequel, we put x = (x,y). 

We suppose that 7^ is a regular triangulation of f2 with triangles r of maximum sidelength h and let 
Sh be a finite-dimensional subspace of C(fi) n H 1 , which, for small enough h and integer r > 2, satisfies 
the approximation property 

inf {\\w - X \\ + h\\w - x\\i} < Ch s \\w\\ s , l<s<r, (3.1) 

xes h 

when w E H s . (In this section C will denote generic constants independent of h.) On the same triangu- 
lation we denote by Sh the subspace of Sh consisting of the elements of Sh that vanish on the boundary 
i.e. Sh = SVi H Hq. Hence, on Sh we have 

inf {|| W - x|| +h\\w- x\\i} < Ch s \\w\\ s , 1 < a < r, (3.2) 

for w E H s n Hq . We will assume that the elements of Sh and Sh are piecewise polynomial functions 
defined on T/,, of degree at most r — 1 on each r E 7h . 

We consider the symmetric bilinear form ao : Hq x Hq — * K defined by 

ao(u, v) := (u, v) + &(Vu, Vu), u,v E Hq , (3-3) 

which is coercive on i?Q x Hq. In addition, we consider the symmetric bilinear form ajv : H 1 x i? 1 — > R 



that is defined by 

a N (u,v) := (u,v) +6(Vu, Vu), u.veH 1 , (3.4) 

and is coercive on iJ 1 x ff. With the aid of an, 8» we define the elliptic projection operators Rh ■ Hq — > 
Sh, Rh ■ H 1 ^ Sh as follows: 

a D (RhW,x) = a D (w,x), Vw e H^, x ^ S h , (3.5a) 
a N {R h w,x) = o,n(w,x), VwG-ff 1 , X&S h - (3.5b) 

As a consequence of (|3.1|> . (|3.2[) and elliptic regularity we have then 

||iy--Rh,u;||fc < C/i s - fe ||tj;|| s , to € iP n # \ 2 < s < r, k = 0, 1, (3.6a) 

and 

\\w - R h w\\ k < C7i s_fe ||u;|| s , weH s ,2<s<r, k = 0, 1. (3.6b) 
We assume that the triangulation 7/j is quasiuniform. Then, the following inverse assumptions hold on 

llxlU < C^-^UII, v x e<? hl (3.7) 
IMIoo^CT*- 1 Hxll, V X eS,, (3.8) 

The same inverse assumptions hold on Sh as well. 

We will also assume that for the elliptic projections we have the following approximation properties 
in the L°° norm: 

Ww-RhwWo, < Cj(h)\\w\\ ri00 , VweW^nH', (3.9a) 
\\w - RhwWoo < C~/(h)\\w\\ ri00 , Vw€W^, (3.9b) 

where j(h) = h r \ log h\ r with f = if r > 2 and f — 1 when r = 2, cf. [S]. Here, (with norm || • ||r,oo) 
denotes the L°°-based Sobolev space on of order r. 
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In the sequel we will use the discrete Laplacian operator Ah ■ H 1 — > Sh, defined for w G H 1 by 

(A h w ) x) = -(Vitf ) Vx), V X ei (3.10) 

By the divergence theorem, it is easy to check that for w G H 2 there holds 

~~ f dw 1 1 ~ 

(A h R h w,x) = (Aw,x)- 7r xds--r(^,x) + -r(RhW,x), X^S h . (3.11) 

Jan on b b 

Hence, if w G -ff 2 with | an = 0, then 

A h R h w = P h Aw+-(R h -P h )w, (3.12) 

where Ph is the L 2 -projection onto Sh- The identity (|3.11|) is proved as follows: For w G H 2 , x ^ Sh, 
using the definition of Rh we have 

(A h R h w, X ) = -{VR h w^x)-\{RhW,x) + \{RhW,x) 
= -(Vu>, V%) - t(w,x) + -r(RhW,x), 



f ^rXds+(Aw,x)-]-(w,x) + ];(RhW,x)- 
Jan on b b 



Denoting the components of v as (u, v), we define the semidiscrete modified Galerkin method as follows. 
We seek r\h ■ [0,T] — > Sh, Uh,Vh : [0,T] — > S/j, approximations to r/,u,v, respectively, such that 

a N (r) ht , <t>) + (Uh x , 4>) + (v hy , <t>) + ((ilhu h )x, 4>) + ((VhVh)y, 4>) = 0, 4> S 
a.D{uh t ,x) + iVh x ,X) + (uhUhx)X) + (vhv hx ,x) + c((&hVh)x,x) = 0, X e 

^ i ^ -/ . (3.13) 

a D (v ht ,i>) + {rihyA) + {uhUhyiii) + (v h v hy ,4>) + c((A h r) h ) y , ip) = 0, ip G Sh, 
r]h(-,0) = RhVo, u h (-,0) = RhU , v h (-,0) = R h v . 

These relations are discrete analogs to the corresponding variational forms of the first p.d.e. of in 
H 1 and of the second in (Hq) 2 . Note that since rjhUh G C(f2), r\hUh\ T G C°°(r) for each r G 7?,,, and 
?7?i 1*^ |an = 0, it follows that rjhUh G Hq. Similarly, r/hVh G -ffp, it 2 G -ffg, w 2 , G Hq. Hence, all terms in 
the inner products of (|3.13p are well defined. We now consider the mappings f x , f y ■ L 2 — > Sh defined for 
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w e L 2 by: 

a>D{fx{w),4>) = (w,<j> x ), <j> e -5ft) 
aD{f v (w),(j>) = (w, (f> y ), <j> 6 5 h , 

and the mappings and g x .g y '■ L 2 — * Sh defined for w G L 2 by 

a;v(<L(w), x) = (w, X e 5ft, 

a.N(g y (w), x) = (w, x y ), x 6 <5/j. 

Then (|3.13p can be written as 

rjh t = F(u h ,v h ,r)h), 

U ht = G(Uh,Vh,Tjh), 

v ht = Z(u h ,v h ,r] h ), 
r) h (0) = R h r]o, Uh(0) 

where 

F(u h ,v h ,r) h ) := g x {u h ) + g y (v h ) + g x (r) h u h ) + g y (i] h v h ), 
G(u h ,v h ,r) h ) := / x (%) + ^ (/*(«;[) + + c fx{^hVh), 
Z(u h ,v h ,rj h ) := fy(r) h ) + ~{f v {u 2 h ) + f y (v 2 h )) + c/ y (A ft %). 

For the mappings / x , f y ,g x and gj, we have the following stability estimates: 
Lemma 3.1 There exists a constant C such that 

(i) ||/xH||i < C\\w\\, and ||/ a H||i < C|M|, u; £ L 2 . 

(ii) II^MIIi < C\\w\\, and \\g v (w)\\i < CHI, w £ L 2 . 

Proof: The proof follows immediately from the coercivity of ao on Hq x Hq and of ajv on if 1 x H 1 and 
the definitions of f x , f y , g x , g y . □ 



< t < T 



(3.14) 



= R h u , v h (0) = R h v a , 
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We define now the negative norms || • ||_i and || • II _2 for functions w G L 2 as 



I II ( W ' Z ) A II II ( W ' 2: ) 

\w\\-i := sup ——r — and M|-2 := sup 



zeH 1 Flu ~ zeH 2 nHi \\ z h 

Lemma 3.2 There exists a constant C > such that 

||/x(x)ll<C||xll-i, ||/„(x)ll <C||xlU v x ei (3.15) 

Proof: Let x 6 Sk- Consider the problem Lw = \x with w — on 90, where L := I — bA with domain 
-ff 2 fl TJq . Then, by elliptic regularity and denoting by L~ 1 the inverse of L, we have 

11 I, (Xx,z) {Lw,z) 

\\Xx\\-2 = SUp I, |, = SUp 



" II II r 11 || 

o^ei? 2 nHj Il z ll2 o^eH 2 ni?J II 2 II 2 

(Lw, HI 2 >r IH 2 

" ||L-%|| 2 ||L-%|| 2 " || W || 

= C\\w\\. (3.16) 



Let now ^ z e H 2 n Hq. Then 

(xx,^) = jx^j < llxll-ill^lli < n I, ^ 



Hence 



Therefore 



and by (gUD 



* a ^2 22 



(Xx,z) 
sup ji < 
OjtzeH 2 nHi \\ z \\2 



||Xx||-2<||xl|-l, 



HI<C|lxll-i- (3-17) 
Consider RhW, the elliptic projection of w onto Sh- Note that Rh,w = f x (x)- In addition, by (|3.6a[) 

\\R h w - w\\ < Ch 2 \\w\\ 2 < Ch 2 \\ Xx \\ < Ch 2 h- 2 \\ X U = C\\x\\-x. 

(In the last inequality we used the inverse inequality ||y||i < C/i~ 2 ||<^||_i, \/tp 6 S/,, which is valid since 
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by (EU): 

IMU= sup ^>ff >c^fjl = ^H|,) 

O^zEffi \\Z\\l \\(p\\l \M\l 

We conclude that 

\\R h w\\-\\w\\<\\R h w-w\\<c\\ x \\-i, 

and by IpTTTji 

||/ x (x)|| = ||JMI < qixll-i, 

which is the required conclusion. The proof for f y is entirely analogous. □ 
Lemma 3.3 There exists a constant C such that 

||A h «||_i < C||Vti||, Vueff 1 . (3.18) 

Proof: Let u G H 1 . Then, for each ic^Oe H 1 we have 

(A h u,w) (A h u,P h w) \\Vu\\\\VP h w\\ ^ 

— n — n — = — n — n — r~n c Vu > 

IfIIi IMIi Hli 

from which, taking the supremum over w S iJ 1 we obtain (|3.18j) . We used the stability of Ph on H 1 , i.e. 
the inequality 

||PH|i <C|M|i, forweH 1 , 

which follows by the argument in [UTj . □ 
We now state and prove the main result of this section. 

Theorem 3.1 For h sufficiently small, the semidiscrete problem \3. 11$ has a unique solution (rjh,Uh,Vh) 
in the interval [0, T] of maximal existence of the solution (r),u,v) of the ibvp HZ\) . Moreover, for some 
constant C — C{n, u, v, T) independent of h we have 

Wv-VhW + \\u-u h \\ + \\v-v h \\ < Ch r , 

and 

\\V ~ Vhh + W u - Ufclli + ll« - "fclli < Ch r -\ 

for each t £ [0, T] . 
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Proof: We suppose that for some constant M there holds that |] 77 1 1 oo < M, ||m||oo < M and ||u||oo < M 
for < f < T. Then, from ([3Jb]) for h sufficiently small 

■fclU < \\Vh - VoWoo + holU = \\RhV0 - m\\oc + IMU < C 1 {h)\\r j0 \\ r , oo + ||ry ||oo < 2M. 



Similar estimates hold for u° and uP. The o.d.e. system ()3.14j) has a unique solution locally in t. By 
continuity, we may assume that there exists th G (0,T] such that Hu/Hoo < 2M, ||^h||oo < 2M and 
IKIIoc < 2M for alH < i h . 
We let now 

p = rj- R h ij, 6 = R h r\ - rj h , t = v - R h v, ( = R h v -v h , a = u - R h u, £ = R h u - u h , 

so that 9 e Sh, (, £ S S/, and 77 — 7/h = p + 0, w — Uh = a + u — w/i = r + By and (|3.14|) we have 

# f = g x {v + + 9y( T + () + 9x{ur/ - u h r) h ) + g y (vr] - v h r/ h ), (3.19) 
6 = f x {e + p) + ^[f x (u 2 )-f x (ul) + Uv 2 )-fM)}+cU^-^hVh), (3.20) 
Ct = / y (r + C) + i{/,(^)-/,(^J + A(^)-/,(«, 2 j}+c/ y (Ar ? -A, l% ). (3.21) 

The equation (|3.20[) may be written as 

= f x (0 + p) + \ {/xKct + 0) + + 0«0 + /aMr + 0) + /x((t + CH)} + c/x(A77 - A h r) h ). 
Taking L -norms and using LemmaEnj $51ty . (|3~T^1 . (|XHa|l . lpT6b"|) and (gH^) we obtain, for < t < t h , 

Mt\\ < \\U8 + p)\\ + \(\\f*M* + Q)\\+ 

\\M<r + t)uh)\\ + ||/x(«(r + C))|| + ||/«((r + CK)ll) + kU&V ~ &hVh)\\ 

< c(\\e + HI + 1|«(<7 + oil + IK* + 0«fcll + Mr + oil + ||(t + CKII + 

\\f x (A v - A h R hV )\\ + \\f x (A h (R hV - 7,011) 

< c[\\e\\ + \\ P \\ + H| L -(H + ||eil) + H|l-(IM| + Hell) + 

IMk-(H + IICII) + Klk~(H + HOI) + IIAr? - A h R hV \\ + \\A h (R hV - %)||-i] 

< C\h r + \\0\\ + HOI + IICII + IIAr? - PftAr?!! + \\R hV - P hV \\ + \\R m - t^Ii! , 
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and so 



Ut\\<c(h r + \\e\\ 1 + u\\ + \\a). 



(3.22) 



Similarly, we have for < t < th 



IMI <^ r + ll^lli + IICII + IICII). 



(3.23) 



The equation (|3.19|) may be written as 



Ot = g x {<J + 0+ 9y(r + C) + 9x(u(p + 8)) + g x ((<j + £)Vh) + g y (v(p + 9)) + g y ((r + ()r, h ). 



Hence, taking 7? 1 -norms and using Lemma l3~Tl and (|3.6b .b) we have for < t < th 



\\9t\h < C-(||a + e|| + ||r + C|| 

+\Hp+9)\\ + \\(<T + OVh\\ + \Hp + 9)\\ + ||(r + C)%||) 

< C(||cr|| + ||e|| + ||r|| + HCll + H|l»(||p|| + \\9\\) + (\\a\\ + M\\)\\Vh\\ L ~ 

+IHk-(IHI + IMI) + (W + ||CII)ll%lk-) 

< c(H + ||a + ||r|| + ||CH + ||p|| + W, 



and therefore that 



^ r + l|0|| + ll£ll + IICII), o<t<t h . 



(3.24) 



From ([3~2"2")) - (|3~2"4"l) we see that for 



< t < t h 



+ U\\ 2 + IICII 2 ) < ch 2r + c(\\e\\t + IICII 2 + IICII 2 ), 



from which, using Gronwall's inequality, we have for < t < th 



ll<?lli + lieii + iicii <chr. 



(3.25) 



Hence, for t < th there holds 



\V - Vh\\ + \\u - u h \\ + \\v- v h \\ < Ch r . 



(3.26) 
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Furthermore for t<t h , we have, by (|X5|) . (f3~9b)) and (jEH} 

\\r]h ~ vWl^ < \\Vh - RhV\\L°° + \\RhV - v\\l°° 

= Ch~ 1 \\9\\ 1 + C 1 {h) 
< Ch r - X 

Therefore, ||%||l«> < ||^||z,«> + ||% — f]\\L°° < CbT^ 1 + M < 2M for h sufficiently small. Similar estimates 
hold for u/j, Vh. These contradict the maximal property of th and we conclude that we may take th = T. 
Hence (|3.26p holds up to t = T, giving the desired optimal-rate L 2 -estimate. The 0(/i r_1 ) H 1 estimate 
follows easily by (j3~6k b) and ([37]). □ 
It is worthwhile to note that temporal derivatives of arbitrary order of the semidiscrete solution 
{rjh,Uh,Vh) are bounded on [0,T] by constants independent of h, as the following proposition shows. 

Proposition 3.1 For h sufficiently small, let (rjh,Uh,Vh) be the solution of the semidiscrete problem 
\3. 14\ ) for t E [0, T] . Then, for j = 0, 1, 2, 3, . . ., there exist constants Cj independent of h, such that 

max (\\di Vh \U + \\d>u h \\ + \\div h \\) < Cj. (3.27) 

Proof: From Theorem 13. II we have, for some constant Co independent of h, 

^t,(Wi+K|| + K||)<C . (3.28) 

Now, from (l3~Ti|) . (ii) of Lemma O and (|3~2^|) . there follows for < t < T 

Wvhth < \\9x(u h )\\i + \\g y (vh)\\i + \\g x (rihUh)\\i + \\g y (vhVh)\\i 

< C(\\u h \\ + \\v h \\ + \\r) h u h \\ + \\vhv h \\) 

<C(1 + ||j?h||oo). (3.29) 
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In addition, from (|3"T4"j) . (i) of Lemma O (|3~T5)) . and (f^5|) we have for < t < T 



\\u ht \\ < \\Uvh)\\ + ^\\Uul)\\ + -\\fM)\\ + \c\\\f x (A hVh )\\ 
<c([W + Kll + Kll + l|A h %||_i) 



< c{\ + IKHooiKii + iK'hllooll^ll + IMIi) 



< C(l + ||ttfc||oo + IK||oc)- 



(3.30) 



A similar inequality holds for ||v/j{||. Now, from the closing argument of the proof of Theorem 3.1 we 
may infer that 



and, consequently, in view of (|3.29|) and (|3.30p . the validity of (|3.27p for j = 1. Differentiating now the 
equations in (|3.14p with respect to t and using again Lemma l3~Tl (13.27P for j = 1, and (|3.3ip we see that 
(pHTf holds for j = 2 as well. 

If we take the second temporal derivative of both sides of the first o.d.e. in (|3.14p . we see that in order 
to obtain a bound for we need, in addition to already established estimates, an h- independent 

bound for ||r7/ llf ||oo- This is obtained as follows: For < t < T we have, for 9 = Rhf] — ijh, from (13. 8p . 
§M, and jS^l, that H^IU < H^IU + \\R h m\\oo < Ch^ptW + C*/(h) + C < C. We similarly get 
^.-independent bounds for ||oo and ||w/it||oo that yield in turn similar bounds for Hc^-u^H and ||<9fv/j||. 
Hence (I3.27P holds for j = 3 too. The case j = 4 follows immediately, as it does not need any L°° bounds 
on temporal derivatives of the semidiscrete approximations of order higher than one. To obtain (|3.27p for 
j = 5 one needs, in addition to already established bounds, ft,- independent bounds for l^, Hc^u/jHoo, 

and || f/i || oo on [0,T]. These may be derived as follows: Differentiate with respect to t the expression 
for 9t after (I3.23P and use the uniform bound on H^Hoo and p.22p ~ (|3.25p to obtain that ||#tt||i < Ch r . 
The required bound for 1 1 c?^ 77^, 1 1 ^ follows then from (|3.8p and (|3.9bp . Similarly differentiating e.g. the 
expression for after (|3.2ip and using the bounds on ||u/i t ||oo) ||t>/ulloo and (|3.22p - (|3.25p we obtain that 
HCttll ^ Ch r , from which H^u/jUqo < C follows. The case j = 6 requires no additional L°° bounds. 

We continue by induction. If j = 2k + 1, L°°-bounds for d^rjh, d^Uh, d^Vh are found by differentiating 
the expressions for 9 t , and £t and using previously established bounds. The even case j = 2k + 2 



max/H^Hoo + \\uh\\oo + IKIU) < C, 



(3.31) 
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requires no additional L°° bounds. (As a corrolary from the above proof it follows also that 

m^llS^lloc + Halloo + Halloo) < C' v 

holds for j = 0, 1, 2, . . ., where Cj are constants independent of h.) □ 
From the result of this proposition we see that the system of o.d.e.'s (|3 . 14(1 is not stiff. Therefore, one 
may use explicit time-stepping schemes to discretize ([3.14)1 in the temporal variable without imposing 
stability mesh conditions on the time step At in terms of h. Error estimates of optimal order in space 
and time in the case of explicit Runge-Kutta full discretizations may be established along the lines of the 
proof of Proposition 10 of |ADMII| . 

Remark 3.1 The H 1 error estimate in Theorem 13.11 may be strengthened as follows. For w G H 1 
define the discrete norm || • \\2.h as 

\\wh, h := (WwWl + WAhwf) 1 / 2 , VweH 1 . 

Then, by letting w £ Hq and considering the boundary- value problem / — bAf = — w x in fl with ^ = 
on dfl, we easily see that g x (w) — Rhf ■ A straightforward computation using (|3.12[l yields then that 
H&cC^OIl! h — Cll^xll 2 ) Il<?2/( w )ll2 h — C\\ w y\\ 2 - We may take now the || • \\2,h norm in (|3.19|) and obtain 

\\e t h,h<c(h r + \\6\\ + u\\ + K\\), o<t<t h , 

instead of (|3.24p . We conclude, along the lines of the proof of Theorem 13. 1[ that 

\\v - Vhh.h + \\u - uhWx + \\v - v h \\i < ChT- 1 . 

4 Numerical experiments 

In this section we present the results of numerical experiments that we performed in the case of the 
Bona-Smith systems using the modified Galerkin method as a base spatial discretization scheme. For the 
temporal discretization of the system of o.d.e.'s (|3.14[) we used an explicit, second-order Runge-Kutta 
method, the so-called "improved Euler" scheme. For the solution of the resulting linear systems at each 
time step we used the Jacobi-Conjugate Gradient method of ITPACK, taking the relative residuals equal 
to 10 for terminating the iterations at each time step. In the computation of the discrete Laplacian 
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A/i we used lumping of the mass matrix. In what follows we present numerical results that confirm the 
expected rates of convergence of the fully discrete scheme and three numerical experiments illustrating 
the use of the method in various surface flows of interest. 

4.1 Experimental rates of convergence 

In order to check the spatial convergence rates of the fully discrete modified Galerkin method we applied 
the scheme to the ibvp (|7^|) for the Bona-Smith system with 2 = 9/11 taking as exact solution 

rj(x, y, t) = cos(7rx) cos(7ry)e*, 
u(x,y,t) = x cos((ttx)/2) sin(vry)e , 
v(x,y,t) — ycos((ny)/2) sin(7rx)e*. 



defined on [0, 1] x [0, 1], which was covered by a uniform mesh consisting of isosceles right-angle triangles 
with perpendicular sides of length h = ^/2/N, where N is the number of the triangles. The solution 
was approximated in space by continuous, piecewise linear functions on this triangulation. We took 
At = 0.01 and computed up to T = 1. In Tables [4~T1 and [4~2l we show the resulting L 2 and H 1 errors and 
the corresponding experimental convergence rates that confirm the result of Theorem 13.11 

Table 4.1: L 2 errors and convergence rates for the modified Galerkin method for the Bona-Smith system 
with 2 = 9/11. Linear elements on triangular mesh and second-order RK time-stepping. 



N 


\\V - Vh\\ 


rate(^) 


\\u-u h \\ 


rate(it) 


\\v - v h \\ 


rate(i>) 


512 


1.3016E-2 




9.0816E-3 




9.0429E-2 




2048 


3.2571E-3 


1.9986 


2.3603E-3 


1.9440 


2.3604E-2 


1.9378 


8192 


8.1485E-4 


1.9990 


5.9882E-4 


1.9788 


5.9882E-4 


1.9788 


32768 


2.0523E-4 


1.9893 


1.5417E-4 


1.9576 


1.5417E-4 


1.9576 



Table 4.2: H 1 errors and convergence rates for the modified Galerkin method for the Bona-Smith system 
with 9 2 = 9/11. Linear elements on triangular mesh and second-order RK time-stepping. 



N 


h - Vh\\i 


rate (77) 


\\u - U h \\l 


rate(it) 


||« - Ufc||l 


rate(w) 


512 


4.3739E-1 




1.8597E-1 




1.8488E-1 




2048 


2.2030E-1 


0.9894 


8.7012E-2 


1.0958 


8.7008E-2 


1.0874 


8192 


1.1038E-1 


0.9970 


4.2013E-2 


1.0504 


4.2013E-2 


1.0503 


32768 


5.5221E-2 


0.9992 


2.0680E-2 


1.0226 


2.0680E-2 


1.0226 
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4.2 Solitary-wave-like pulse hitting a cylindrical obstacle 

In our first experiment the domain that we consider is the rectangular channel [—15, 15] x [—30, 50] in 
which is placed a vertical impenetrable cylinder centered at (0,10) with radius equal to 1.5. We pose 
normal reflective boundary conditions on the boundary of the cylinder and along the lines y = — 30 and 
y = 50, and homogeneous Neumann boundary conditions for rj and v on the lateral boundaries x = ±15. 
As initial conditions we use the functions 



T] (x,y) = Asech 




u (x,y) = 0, (4.32) 
vo(x,y) = r)o(%,y) - \ril[x,y), 

where A = 0.2 and c s = 1 + 4, which represent a good approximation to a line solitary wave, (DMSIj . 
of the BBM-BBM system. We integrate under these initial and boundary conditions the Bona-Smith 
system (|l.lj) with 8 2 = 9/11 and compare its evolving solution with that of the BBM-BBM system, i.e. 
(fi~T|l with 9 2 = 2/3. The Bona-Smith system was discretized in space by the modified Galerkin method 
(amended in a straightforward way to handle the Neumann conditions on v on the lateral boundaries) 
with continuous, piecewise linear elements on a triangulation consisting of N = 290112 trangles; this 
mesh was fine enough for the 'convergence' of the numerical solution. The same truangulation was used 
to discretize in space the BBM-BBM system with the standard Galerkin method, DMSI . Both systems 
were discretized in time with the improved Euler method with At = 0.01. 

The line solitary-wave-like pulse is centered at y = —10 at t = 0. It propagates, mainly in the positive 
y direction, travelling with a speed c s of about 1.05, cf. Figure 4.2. The bulk of the solitary wave travels 
past the cylinder, while smaller amplitude scattered waves are produced by the interaction of the wave 
with the obstacle. Figure 4.3 shows, for both systems, contours of the elevation of the wave, superimposed 
on velocity vector plots, near the obstacle during the passage of the solitary wave. Figure 4.4 shows the 
free surface elevation for both systems as a function of y at x = at several temporal instances. It 
is worthwhile to note that the maximum run-up at the extreme upstream point (x,y) = (0,8.5) of the 
cylinder was equal to z = 0.335107 (achieved at t — 16.72) for the BBM-BBM system; the corresponding 
run-up for the Bona-Smith system was z = 0.319960 at t — 16.74. The analogous values for the run-up on 
the downstream point (0, 11.5) of the cylinder were z = 0.234310, t = 22.71 for the BBM-BBM system, 
and z = 0.230590, t = 22.68 for the Bona-Smith system. Figure 4.5 shows the history of the free surface 
elevation for both systems at the extreme points y — 8.5, y = 11.5 along the x = diameter of the 
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cylinder, perpendicular to the impinging wave. In general, we did not observe great differences in the 
behaviour of the solutions of the two systems. 

4.3 Evolution and reflection at the boundaries of a 'heap' of water 

The sequence of plots of Figure 4.6 shows the temporal evolution of the free surface elevation of an initial 
Gaussian 'heap' of water with 

V0 (x,y) = 2e -^+ 4 °) 2 +fe+ 40 > 2 )/ 5 , y (x,y) = 0, 

as it collapses, forms radial outgoing riples and interacts with the reflective boundaries near a corner of 
the square [—80, 80] x [—80, 80]. Again, no major differences were observed between the solutions of the 
two systems (BBM-BBM and Bona-Smith, 6 2 = 9/11). The computation was effected with N = 84992 
triangular elements and At — 0.1. Figure 4.7 shows the corresponding one-dimensional plots (along 
x = —40) of the free surface elevation for both systems at four temporal instanses. 

4.4 Line wave impinging on a 'port' structure 

In this experiment we integrated the BBM-BBM and the Bona-Smith (9 2 = 9/11) systems in dimen- 
sional variables. Our domain represents part of a 'port' of depth h — 50m consisting of the rectangle 
[-250,250] x [0,2000] minus the rectangular 'pier' [0,100] x [0,700]. (All distances in meters). Normal 
reflective boundary conditions are assumed to hold along the boundary of the pier and the intervals 
[—250, 0] and [100, 250] of the x-axis, while homogeneous Neumann boundary conditions have been im- 
posed on rj and v on the remaining parts of the boundary. An initial wave with rj (x,y) = Ae~ y sooo° , 
A = lm, and Uo(x,y) = Orn/sec, v (x,y) = —\{r)o{x,y) — jTIq(x, y)) m/sec travels mainly towards the 
negative y direction shedding a dispersive tail behind. (Both systems were integrated with N = 77632 
triangles and At = 0.01 sec). For the impinging wave at t = 30 sec we estimated that A/ho — 0.014, 
X/ho = 16.6, so that the Stokes number is approximately equal to 3.9, within the range of validity of 
the Boussinesq systems. The incoming wave hits the pier and the port boundary at x — 0, is reflected 
backwards and interacts with the remaining boundary. Figure 4.8 shows contour plots of the free surface 
elevation for both systems at several temporal instances in the whole domain. In Figure 4.9 we plot the 
computed free surface elevation as a function of y along the x — 40 m line at several temporal instances 
as the wave hits the pier front (y = 700) and reflects backwards. Most of the differences in the solution of 
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BBM-BBM t = 40 



Figure 4.2: Experiment 4.2 Free surface elevation 
{6 2 = 9/11) systems. 




Bona-Smith t = 10 




Bona-Smith t = 20 
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four time instances. BBM-BBM vs. Bona-Smith 
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t = 30 t = 40 

Figure 4.4: Experiment 4.2. Free surface elevation plots as functions of y at x = at five time instances. 
BBM-BBM , Bona-Smith — . 
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the two systems are of the order of 10 cm and are observed in the reflection phase. Figure 4.10 shows the 
temporal history of the free surface elevation for both systems at the point (a;, y) = (43.75, 700) at the 
front of the pier. The maximum run-up observed at that point was equal to z = 0.976 m (at t = 38.6 sec) 
for the BBM-BBM system and to z = 0.988 m (at t = 38.7) for the Bona-Smith system. 
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BBM-BBM t = 80 Bona-Smith t = 80 

Figure 4.8: Experiment 4.4. Free surface elevation at four time instances. BBM-BBM and Bona-Smith 
(8 2 = 9/11) systems, (elevation and x,y in meters). 
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Figure 4.9: Experiment 4.4. Free surface elevation (in meters) as function of y at four time instances 
along x = 40m. BBM-BBM , Bona-Smith (6» 2 = 9/11) — . 
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